#include <iostream>
#include <string>
#include <vector>
#include "getopt_pp.h"

#define WGDALRASTER_H_IMPLEMENTATION 1
#include "wGdalRaster.h"

#define W_GEOJSON_AND_SHP_READER_IMPLEMENT 1
#include "wgeojsonandshpreader.h"

#define W_LEVEL_HSEG_IMPLEMENT 1
#include "wlevelhseg.h"

using namespace std;
using namespace GetOpt;

int main(int argc,char* argv[]) {

	cout<<"A program to make mask raster file from lon lat and vector files."<<endl;
	cout<<"v1.0.1.2d 2023-11-10"<<endl ;
	cout<<"usage:"<<endl;
	cout<<"vector2orbitmask --vec vectorfile --reso resolution [--lon lonfile --lat latfile]|[--m1b modisl1bfile]|[--m03 modisgeofile]|[--fgeo fygeofile] -o maskout.tif [--size xsize ysize]"<<endl;
	cout<<" --vec geojson or shp polygon, should be wgs84 crs."<<endl;
	cout<<" --reso resolution in wgs84"<<endl;
	cout<<" --lon longitude file, 1st order"<<endl;
	cout<<" --lat latitude file , 1st order"<<endl;
	cout<<" --m1b modis l1b file, 2nd order"<<endl;
	cout<<" --m03 modis geo file, 3rd order"<<endl;
	cout<<" --fgeo fy3d geo file, 4th order"<<endl;
	cout<<" --size new output size, optional"<<endl;
	cout<<" -o output mask file 1=in, 0=out"<<endl;
	cout<<"-------"<<endl;
	GDALAllRegister();

	vector<int> newSize ; 
    string vectorfile ;
    double resolution=0;
    string lonfile ;
    string latfile ;
    string outfile;
    string m1bfile;
    string m03file;
    string fygeofile;

    GetOpt_pp ops(argc,argv) ;
    ops>>Option("vec",vectorfile);
    ops>>Option("reso",resolution);
    ops>>Option("lon",lonfile) ;
    ops>>Option("lat",latfile) ;
    ops>>Option('o',outfile) ;
    ops>>Option("size",newSize) ;

    ops >> Option("m1b" , m1bfile ) ;
    ops >> Option("m03" , m03file ) ;
    ops >> Option("fgeo" , fygeofile ) ;

    cout<<"vec "<< vectorfile<<endl;
    cout<<"reso "<< resolution<<endl;
    cout<<"lon "<< lonfile<<endl;
    cout<<"lat "<< latfile<<endl;
    cout<<"o "<< outfile<<endl;
    if( newSize.size()==0 ) {
    	cout<<"size "<< "none"<<endl;
    }else if( newSize.size()==2 ) {
    	cout<<"size "<< newSize[0]<<" "<<newSize[1] <<endl;
    }else{
    	cout<<"size should be two int values, or not use."<<endl;
    	return 12 ;
    }
    cout<<"m1b "<< m1bfile<<endl;
    cout<<"m03 "<< m03file<<endl;
    cout<<"fgeo "<< fygeofile<<endl;
	cout<<"-------"<<endl;

	if( vectorfile=="" ) {
		cout<<"vectorfile is empty"<<endl;
		return 11 ;
	}

	if( resolution < 0.000001 ) {
		cout<<"invalid resolution"<<endl;
		return 12 ;
	}

    if( lonfile=="" && latfile=="" ) {
    	if( m1bfile!="" ) {
    		lonfile = string("HDF4_SDS:UNKNOWN:")+m1bfile+":1" ;
    		latfile = string("HDF4_SDS:UNKNOWN:")+m1bfile+":0" ;
    	}else if( m03file!="" ){
			lonfile = string("HDF4_SDS:UNKNOWN:")+m03file+":1" ;
    		latfile = string("HDF4_SDS:UNKNOWN:")+m03file+":0" ;
    	}else if( fygeofile!="" ) {
    		lonfile = string("HDF5:")+fygeofile+"://Geolocation/Longitude" ;
    		latfile = string("HDF5:")+fygeofile+"://Geolocation/Latitude" ;
    	}else {
    		cout<<"no lon,lat,m1b,m03,fgeo"<<endl;
    		return 12 ;
    	}
    }

    wGdalRasterFloat lonRas,latRas;
    bool lonok = lonRas.open(lonfile);
    bool latok = latRas.open(latfile) ;
    if( lonok==false ) {
    	cout<<"failed to open lonfile "<<lonfile<<endl ;
    	return 13 ;
    }
    if( latok==false ) {
    	cout<<"failed to open latfile "<<latfile<<endl ;
    	return 13 ;
    }
    int llxsize = lonRas.getXSize() ;
    int llysize = lonRas.getYSize() ;
    if( llxsize!=latRas.getXSize() || llysize!=latRas.getYSize() ) {
    	cout<<"lon lat size not same."<<endl;
    	return 14 ;
    }

    vector<WLineSeg> lineSegs; 
    string error;
    double vecx0,vecx1,vecy0,vecy1;
    bool vecok = wGeoJsonShpReader::readFromFile2(vectorfile,lineSegs,vecy0,vecy1,vecx0,vecx1,error) ;
    if( vecok==false ) {
    	cout<<"read vector file failed: "<<vectorfile<<endl;
    	cout<<error<<endl;
    	return 15 ;
    }

    cout<<"vector bounding(x0 x1 y0 y1):"<<vecx0<<" "<<vecx1<<" "<<vecy0<<" "<<vecy1<<endl;
    cout<<"lineSegs "<<lineSegs.size()<<endl;
    if( lineSegs.size()==0 ) {
    	cout<<"lineSegs is empty."<<endl;
    	return 16;
    }

    int fakeWidth = (vecx1 - vecx0) / resolution ;
    if( fakeWidth==0 ) {
    	cout<<"vec x0 x1 lower than resolution can not calculate hseg "<<vecx0<<" "<<vecx1<<" "<<resolution<<endl;
    	return 17;
    }
    wLevelHseg hseg ;
    bool hsegok =  hseg.init2(
        lineSegs, //矢量文件所有线段
        vecy0 ,
        vecy1 ,
        vecy0 ,
        vecy1 ,
        vecx0 ,
        fakeWidth ,
        resolution, //栅格分辨率，使用正值
        error
        ) ;
    if( hsegok==false ) {
    	cout<<"failed to compute hseg:"<<error<<endl;
    	return 18 ;
    }
    cout<<"hseg count : "<<hseg.hsegs.size()<<endl;

    vector<unsigned char> maskdata(llxsize*llysize,0) ;
    int goodcnt = 0 ;
    for(int iy=0;iy<llysize;++iy) {
    	for(int ix = 0 ; ix < llxsize; ++ ix ) {
    		float lonval = lonRas.getValuef(ix,iy,0) ;
    		float latval = latRas.getValuef(ix,iy,0) ;

    		int hx = (lonval - vecx0)/resolution;
    		int hy = (vecy1 - latval)/resolution;

    		int findy = hseg.binarySearhHsegIndexByY(hy) ;
    		if( findy<0 ) continue;
    		for(int hy2 =findy; hy2 < hseg.hsegs.size();++hy2 ) {
    			if( hy2!=hy )break;
    			if( hx>= hseg.hsegs[hy2].x0 && hx <= hseg.hsegs[hy2].x1 ) {
    				++goodcnt;
    				maskdata[iy*llxsize+ix] = 1 ;
    				break;
    			}
    		}

    	}
    }
    cout<<"good count: "<<goodcnt<<endl;
    if( newSize.size()==0 ) {
    	newSize.resize(2);
    	newSize[0] = llxsize;
    	newSize[1] = llysize;
    }
	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* dsout = driver->Create( outfile.c_str() , newSize[0]  , newSize[1]  ,1,GDT_Byte,0) ;
	dsout->GetRasterBand(1)->RasterIO(GF_Write,0,0,newSize[0],newSize[1],maskdata.data() ,llxsize,llysize,GDT_Byte,0,0,0) ;
	GDALClose(dsout) ;

	cout<<"done"<<endl;
	return 0 ;
}